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Abstract 

We examine the two-dimensional extension of the model of Kesslcr and Sander of competition between 
two species identical except for dispersion rates. In this class of models, the spatial inhomogeneity of 
reproduction rates gives rise to an implicit cost of dispersal, due to the tendency to leave favorable 
locations. Then, as in the Hamilton-May model with its explicit dispersal cost, the tradeoff between 
dispersal case and the beneficial role of dispersal in limiting fluctuations, leads to an advantage of one 
dispersal rate over another, and the eventual extinction of the disadvantaged species. In two dimensions 
we find that while the competition leads to the elimination of one species at high and low population 
density, at intermediate densities the two species can coexist essentially indefinitely. This is a new 
phenomenon not present in either the one-dimensional form of the Kcssler-Sander model nor in the 
totally connected Hamilton-May model, and points to the importance of geometry in the question of 
dispersal. 

Introduction 

A central question in ecology is how species determine their migration rate. Migration entails a number 
of possible benefits to the population, among them reducing overcrowding areas and the discovery of 
new resources to exploit. There are also significant costs associated with migration. Especially in patchy 
environments, there may be a significant risk of death (due to predation or other exogenous causes) in 
trying to cross over unsuitable areas to find a colonizable site. Presumably, migration rates are under 
evolutionary control, and so can in principle be adjusted to some optimal level. There have been a large 
number of studies in this field, both from an analytical and a numerical perspective. 
The current work is an extension of the work of Kcssler and Sander [1] reconciling two seemingly contra- 
dictory approaches to this question of optimal dispersal. On the one hand, there has been extensive study 
of the rate equations approach to modeling. This approach had its beginnings in a result of Hastings [2] , 
who showed both in a two patch system and in a continuous one-dimensional system that the species 
with the slower dispersal rate is stable with respect to the introduction of a faster dispersing species, 
always driving the faster species to extinction, as long as the two patches were not identical. This is 
due to the fact that the overall effect of dispersal was to take individuals from favorable areas to less 
suitable places, leaving the slow species in control of the favorable site. The higher production rate of 
progeny at the favorable site eventually leads to a takeover of the poorer site as well. Thus, dispersal in a 
inhomogeneous environment carries with it an implicit cost. This finding was reinforced by Dockery and 
coworkers [3]. Examining the continuous one-dimensional version of the Hastings model, they studied 
the reaction-diffusion system: 

f = Dfr + a{x)f-f3f{f + s) 

s = Dss" + a{x)s ~ Ps{f + s) ^ ' 

Here, / and s represent the local population density of the fast (/) and slow (s) migrating species, with 
Df > Ds ■ The inhomogeneous distribution of resources is modeled by a spatially varying growth rate, 
a(a;), common to the two species. The competition is of logistic type, with each species competing equally 
against both fellow members of the same species and members of the other variety. The Hastings result 
is that the pure slow state with / = is linearly stable, and the pure fast state, s = is unstable. 
Dockery, et al. were able to prove the stronger result that, independent of the exact form of a(x) and 
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of the initial conditions, the pure slow state is the only global attractor and the faster species always 
went to extinction. The implication of both these works is that any mutation resulting in a lower level 
of migration should fix in the population, driving the migration rate eventually to zero. 

This strong result is in sharp contrast to what occurs in the individual based models originated by 
Hamilton and May |4]. Looking at a system with an infinite site of patches fully connected by dispersal, 
and with an explicit cost of dispersal, they showed that the evolutionary stable stategy was a large level 
of dispersal. This optimal level of dispersal approached unity as the dispersal cost was lowered to zero, 
and surprisingly, remained finite even when the dispersal cost approached 100% of the dispersers. 

Kessler and Sander [T] reconciled these results by showing that the demographic fluctuations inherent 
in the individual based model, and absent in the rate equation approach, were responsible for the very 
different results of the two models. By adding a noise term to the differential equation. Eq. ([T]), they 
showed that fluctuations have a negative effect and that dispersal, by reducing such fluctuations by 
uniformizing the population in space, has a beneficial effect. Thus, the tradeoff between the dispersal 
cost and the fiuctuation reduction benefit gives rises to a nonzero optimal value of dispersal rate. This 
result underlies the decrease of the optimal dispersal rate with increasing population density seen in 
the Hamilton-May-Comins [5] extension of their model to arbitrary population density, since increasing 
population density decreases fluctuations and so reduces the beneficial aspects of dispersal. 

To verify their findings, Kessler and Sander studied an individual based version of the one-dimensional 
Dockery model. The system has two species of individuals, fast and slow, residing on sites of a one- 
dimensional lattice = 0, 1, . . . i with no-fiux boundary conditions. The system is designed to reproduce 
a spatially discretized version of the Dockery model in the deterministic limit. All individuals have a 
probability of a{xi)dt of reproducing in some small time interval dt, depending on their location Xi. 
Furthermore, individuals at Xi have a concentration-dependent probability of death, (/^ -I- Si)l3dt, where 
fi and Si are the numbers of fast and slow individuals sharing the particular site. They also move to a 
nearest- neighbor site with probability Df^^dt, depending on their identity. Kessler and Sander chose to 
examine a distribution of resources corresponding to two "oases" of food separated by a desert: 



They set /3 = — , so that in the absence of diffusion, the total population at each site would saturate 



at the local carrying capacity Noa{xi). The results of this model showed that, for a given level of 
inhomogeneity rj, when Nq was large, the rate equations correctly predict that the slower species wins. 
On the other hand, for TVq small, the ability of the faster species to use diffusion to minimize fluctuations 
led to its victory, against the predictions of the rate equations. One surprise was that the boundary 
between large and small Nq was surprisingly large, in the range of hundreds per lattice site, against 
the naive expectation that the rate equations would be reliable for A^o's of order 10 or so. The interplay 
between inhomogeneity and Nq was consistent with the observation that in the Hamilton-May-Comins [5] 
extension of the model, wherein an arbitrary carrying capacity is allowed, the evolutionary stable dispersal 
rate decreases with A^o, and with increasing dispersal cost. 

The model we examine herein is the natural two-dimensional generalization of the Kessler-Sander 
model. The system consists of a square two-dimensional lattice, Xi^yi = 0,1,... L. The dynamics is 
exactly the same as in the one dimensional case, with the dispersal move now occurring to one of the 
four nearest-neighbor sites. We choose for the inhomogeneity a "checkerboard" pattern 




(2) 



1 



No 




(3) 



with oases in the corners and the center and deserts in-between, see Fig. [T] In all our simulations, we 
have set L = 50. 
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Results 




Figure 1. Plot of the local birthrate a{x, y) in our 50 x 50 system. 

Figs. [2]|4] display some snapshots from the two dimensional simulation. The color represents the 
relative abundance of fast and slow on a site, with all fast being pure red and all slow being pure blue. 
Figl2] displays a system with the parameters: ry = 0.4, Nq = 400. Fig. [2]^a) displays the initial state of 
the system for all the runs, namely only slows in the left half of the system and only fasts in the right 
half. In FiglUb)-FigI2t^e) it is possible to see that initially the fasts predominate in the deserts and the 
slows in the oases. Eventually, for this set of parameters, with its relatively low population density, the 
fast dispersing species encroach on the oases and drive the slows toward extinction, as seen in n FiglJl^f). 

Fig. [3] displays a system with the same degree of inhomogeneity, rj = 0.4 but a larger population 
density TVq = 700. The system starts as in FiglJJa). Again, the fasts fare better in the deserts, but the 
eventual advantage accrues to the slows, as expected due to the higher population density. By Figl^Jd), 
after 4500 steps, the victory of the slows is essentially complete, with the slows in the "deserts" making 
a valiant but futile last stand, in accord with the statement that the slows win at sufficiently large A^o- 

The above results are consistent with the picture that emerged from Kessler and Sander's simulations 
of the one-dimensional model. However, we show in Fig. |4] snapshots of a system with the same inhomo- 
geneity rj = 0.4 and an intermediate value of the population density, A^o = 600. Again, the system starts 
as in Fig|2Ja). In this figure it is clearly seen that despite the much longer run times compared with the 
above two cases, there is essentially no change in the ratio between the fasts and slows after 15000 time 



Figure 2. Snapshots from a simulation with the parameters 77 = 0.4, A^o — 400. (a) The initial state of 
the system, (b)-(e) The fasts initially fare better in the "desert" and the slows in the "oases", (b) the 
system after 100 steps, (c) 450 steps, (d) 1000 steps, (e) 1300 steps, (f) Essentially total victory of the 
fasts after 3000 steps. 



Figure 3. Snapshots from a simulation with the parameters 77 = 0.4, A^o — 700. (a) The system after 
100 steps, (b) 200 steps, (c) 1000 steps, (f) almost complete victory of the slows after 4500 steps. 
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Figure 4. Graphic simulation witli the parameters i] = 0.4, A^o = 600. (a) the system after 1000 steps, 
(b) the system after 5800 steps, (c) the system after 11800 steps, (f) the system after 15000 steps, there 
are no essential changes during the time. 

steps. This behavior leads to the speculation that the state displayed in this figure will be preserved also 
at infinite time, with the system having settled into a state of stable coexistence. 

FigIS] display the percentage of fast individuals in the entire system as a function of time for the three 
cases presented in the snapshots. For rj = 0.4, = 400 the percentage of fasts rises quickly to 100%. 
While rj = 0.4, Nq = 700 the percentage of fasts decreases to zero. It should be remarked that in this 
case the victory of the slows in clear from the spatial snapshot long before the total number of fasts has 
dropped significantly. For rj = 0.4, Nq = 600 the percentage of fasts fluctuates around 67.2%, with no 
discernible secular trend. 

To investigate further, we scan the rj, Nq parameter space, the results of which are summarized in 
FigEl The blue color represents parameters in which the slows win, the red color represents parameters 
in which the fasts win. The green color represents parameters with apparent coexistence. The overall 
trend of the dependence on rj and Nq are as in the one dimensional model with higher inhomogeneity 
and higher Nq both favoring the slows. The new feature is of course the intermediate coexistence phase, 
which is seen to shrink in extent as rj increases, at least for rj > 0.4. 

In order to clarify if the apparent coexistence phase is real or only the result of the finite running time 
of the simulation, we measured the time until victory as the coexistence region is approached from either 
above (high Nq) or below. Figs. [7] and |S] display the average time until one species wins as a function 
of Nq. Fig. [7] displays the average time until the fasts win in the small Nq regime and Fig. |8] displays 
the average time until the slows win in the higher Nq regime. The time till extinction of the loser grows 
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r|=0.4 




Figure 5. Percentage of fast individuals as a function of time for the simulations presented in Figs. [21 
El and a with r] = 0.4 and iVo = 400, 600, 700, respectively. 

markedly as A'o approaches the coexistence area, and is consistent with a power-law divergence: 

tf = Bf{Nf,,-No)-^^ (4) 

Likewise, the time to extinction of the losing fast species grows as Nq approaches the coexistence area 
from above, again consistent with a power-law divergence: 

= B,(iV,-7V,,e)-^= (5) 

The parameters of Eqs. (jH O were found by least-squares fitting. The graphs lead to the conclusion that 
also for infinite time there will be no victory of one population because the winning times diverge as the 
coexistence region boundary is crossed. 

It is interesting to note that the exponents characterizing the divergence are different on the two sides 
of the transition. In addition, the exponents decrease with increasing rj. Graphs of the exponents Af, As 
as a function of ry are presented in Fig. [S] 

To point out the contrast between the coexistence dynamics we see in the two dimensional system 
and the behavior of the one dimensional simulation, we have measured the average extinction time as 
a function of A'o for a fixed 77 = 0.6 in the one dimensional system. The data is presented in Fig. [TOl 
Here the extinction time rises as the border between the fast dominated and the slow dominated phases 
is approached, but it shows no sign of diverging. 

Up till now, we have focused only on one specific pattern of inhomogeneity in the two dimensional 
system, namely the " checkerboard" . We now investigate a "striped" pattern of inhomogeneity, which is 
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Figure 6. Results of the two dimension simulation. Red: fasts win, blue: slows win, green: coexistence 



astripo(a:,y) = 1 + J/cos (6) 



a uniform extension in the second dimension of the one dimensional a{x) (Fig. Ilip 

( 

We present snapshots from the simulation in Fig. [12] for ij — 0.4, Nq = 600, which showed coexistence 
in the checkerboard geometry. As always, initially the fasts dominate the "deserts" and the slows the 
"oases" , but as in the checkerboard geometry, the system settles into a coexistence state. There is still a 
correlation between the local growth rate and the fraction of fasts in this coexistence state, but as we see 
in Fig. 1131 it is surprisingly weak, with only a 15% variation between the fast/slow balance between 
the deserts and the oases. 

Given that coexistence behavior is not present in the one-dimensional, i.e. 50 x 1, system, but does 
obtain in the two-dimensional, i.e. 50 x 50, system, the critical parameter underlying coexistence must 
be the width. We show in Fig. [T3]the average time to victory for systems of width 1, 10 and 50. We 
see that as the width increases, the peak moves toward lower Nq and increases in height. The maximum 
time to victory for the 50 x 10 system is 29000 steps, which is much larger than the 50 x 1 system, but 
nevertheless is finite. We conjecture that this should be the case for any finite width system, with the 
maximum time to victory diverging with the system size. For the 50 x 50 system, such a maximum time 
is evidently much larger than we can afford to measure. 



Summary and Discussion 

In conclusion, we have shown that a checkerboard pattern of inhomogeneity leads to the coexistence of 
two different dispersal strategies, for some intermediate range of densities. This feature is absent in both 
one-dimensional systems and in fully-connected systems. The transition to coexistence is marked by the 
power-law divergence of the extinction time of the under-performing strategy. More study is required to 
clarify the nature of this coexistence phase, and the critical behavior associated with the transition. 




Figure 7. The time until the extinction of the slows leaving the fasts victorious as a function of Nq. 
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Figure 8. The time until the extinction of the fasts leaving the slows victorious as a function of Nq. 
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Figure 9. The powers Af, As of the fit to the functional form Eq. (|4l5p . (a) Fasts win, (b) Slows win. 

The phenomenon of coexistence in an inhomogeneous environment is not surprising from the perspec- 
tive of niche theory. The "oases" , with their high reproduction rate, arc evidently favorable for the slows 
and the "deserts" for the fast dispersers. However, given the density-dependent nature of the effective fit- 
ness, this is a simplistic view. In one dimension, one species may exhibit a temporary advantage, however 
in the fullness of time, one species always drives its competitor to extinction. Evidently, however, in the 
two dimensional environment, this is no longer true near the region where the advantage changes between 
species in a one dimensional patterned environment. The question of coexistence evidently requires a 
more refined understanding. 

Even with a standard density-independent fitness landscape, the interplay of dispersal and niches 
raises interesting questions worthy of further study. Dispersal will at a minimum tend to wash out 
the niche boundaries, with some degree of mixing at the niche boundaries. The question that arises is 
whether there is some degree of dispersal that wipes out completely the underlying niche structure of the 
environment and leads to domination by a single "compromise" strategy. In our system, the fairly small 
difference between the population balance in the "oases" and the "deserts" means that dispersal effects are 
strong enough to wash out most, but not all, of the differences between the niches. Given the absence of 
a clear boundary between the two niches, it is surprising that in the two dimensional system, coexistence 
is still possible. Comparing the present model to one with density-independent fitness differences should 
be very instructive. 
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Figure 10. The average time until one of the population wins in the one-dimcnsional system, with 
T] — 0.6. 
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Figure 11. Plot of the local birthrate a{x, y) in the striped 50 x 50 systeni. 
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Figure 12. Graphic simulation of the striped pattern with the parameters 77 = 0.4, Nq = 600. (a) The 
initial state of the system, (b) the system after 200 steps, (c) the system after 400 steps, (d) the system 
after 1000 steps, (e) the system after 5000 steps, (f) the system after 20000 steps, in the (d)-(f) 
snapshots there are no essential changes during the time. 
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Stripes r|=0.4 N=600 




Figure 13. Percentange of fast individuals as a function of position in the striped two-dimensional 
system. 77 = 0.4, No = 600, t = 20000. 
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Figure 14. The average time until one of the population wins in the striped 1 x 50,10 x 50 and 50 x 50 
system, with 77 = 0.6. 



